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In  this  report,  we  discuss  the  development  of  a  multi-scale  framework  to  simulate  and  predict  shear  initiated  reactions  in 
energetic  and  reactive  materials.  First,  we  implemented  the  framework  into  an  Eulerian  wave  propagation  code.  Then,  using  the 
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1.  Objective 


The  objectives  are  as  follows: 

1.  Formulate  and  develop  a  multi-scale  approach  for  simulating  shear  initiated  reactions  that 
will  span  from  the  mesoscale  to  the  full  continuum  scale,  to  link  the  length  scale  of 
material  heterogeneities  to  the  length  scale  of  the  weapon  system. 

3.  Develop  models  and  establish  bridging  at  the  various  length  scales,  followed  by 
verification  and  assessment  of  the  multi-scale  approach. 

The  final  product  is  an  improved  predictive  capability  implemented  into  the  Combined  Hydro 
and  Radiation  Transport  Diffusion  (CTH)  hydrocode  developed  by  Sandia  National  Laboratories 
that  allows  for  the  prediction  and/or  simulation  of  high  explosive  (HE)  and  reactive  material 
(RM)  reactions  when  subjected  to  loads  that  result  in  shear  localizations. 


2.  Approach 


A  schematic  of  the  approach  developed  is  shown  in  figure  1.  It  consists  of  four  major  steps  as 
follows: 

1 .  Perform  the  mesoscale  modeling. 

2.  Convert  the  mesoscale  output  to  input  for  the  reactive  bum  model. 

3 .  Link  the  reactive  bum  model  output  to  the  continuum  level  simulation  of  shear  initiated 
reactions. 

4.  Assess  the  overall  approach  by  implementing  it  into  the  continuum  mechanics  code  CTH, 
and  perform  simulations  for  HEs  and  RMs. 
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As  depicted  in  figure  1(a),  each  mesoparticle  represents  several  atomic-level  unit  cells  of  the 
crystalline  material.  The  mesoscale  method,  Dissipative  Particle  Dynamics  with  Energy 
Conservation  (DPDE)  ( 1 ,  2),  provides  the  input  to  the  reactive  bum  model.  Thermodynamic 
properties  (pressure,  density,  temperature,  and  energy)  from  the  mesoscale  simulations  are 
generated  at  various  steady-state  shear  rates,  providing  look-up  data  tables  for  use  in  the  CTH 
modeling  (cfi,  figure  1(b)).  For  each  strain  observed  in  the  computational  domain  where  a  shear 
localization  has  developed,  the  state  variables  and  the  energy  spectrum  obtained  by  the 
mesoscale  model  are  known  at  the  continuum  level.  With  this  approach,  the  assumptions  made  in 
the  continuum  mechanics  code  are  circumvented.  In  other  words,  by  incorporating  mesoscale 
level  calculations  into  the  continuum  mechanics  code,  more  accurate  estimates  for  the 
thermodynamic  state  within  the  localized  regions,  and  consequently,  more  accurate  input  to  the 
reactive  bum  model  are  obtained. 

We  developed  a  sub-grid  model  to  account  for  energy  dissipation,  which  brings  the  mesoscale 
temperature  estimates  to  the  continuum  level  in  a  relatively  cell-size  independent  manner  by 
calculating  a  volumetric  averaged  solution  of  the  extent  of  reaction  at  the  sub-grid  level  to 
represent  the  continuum  extent  of  reaction  (cf.  figure  1(c)).  We  use  this  extent  of  reaction  in  the 
equation  of  state  to  determine  the  amount  of  energy  release. 

The  sub-grid  calculations  are  performed  for  cells  that  are  identified  to  contain  shear  bands  at  the 
continuum  level.  A  numerical  framework  for  nucleating  and  propagating  shear  band  localizations 
already  exists  in  CTH  (3-5).  When  a  set  of  nucleation  criteria  is  satisfied,  a  shear  band  is  fonned 
by  introducing  Lagrangian  particles,  which  conform  to  local  planes  of  maximum  shear  as  they 
propagate  in  three-dimensional  space  until  a  set  of  growth  criteria  are  no  longer  satisfied  along 
the  points  defining  its  boundary.  We  use  this  framework  to  track  shear  localizations,  and  we 
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applied  our  multi-scale  approach  at  these  locations  to  simulate  the  shear  initiation  of  energetic 
and  reactive  materials  (cf.  figure  1(d)).  We  then  implement  the  approach  into  CTH  and  conduct 
simulations  for  TNT  to  verify  the  algorithm.  The  details  of  the  approach  follow. 

The  energy  conserving  version  of  the  dissipative  particle  dynamics  method  (DPDE)  is  a  particle- 
based  mesoscale  method  that  simulates  the  hydrodynamic  behavior  of  materials,  conserving  both 
momentum  and  energy  while  allowing  the  mesoparticles  to  exchange  both  viscous  and  thennal 
energy  (1,2).  In  the  DPDE  method,  the  changes  with  respect  to  time  t  in  position  {r,}, 
momentum  {p,  },  and  internal  energy  { u,  }  of  the  /th  particle  with  mass  m,  due  to  the  interaction 
with  the  /th  particle  can  be  expressed  as 


P, 


in , 


pl=F#c+F|“‘+F?+F; 


U:  ~  - 


+F,;)+<fy  +q‘ 


(1) 

(2) 
(3) 


with  the  requirement  that  q'.  =  -qi  and q*  =  -qRn ,  where  q'-  and  qR  are  the  dissipative  and 

random  mesoscopic  heat  flows,  respectively,  and  where  the  dot  notation  denotes  time- 
differentiation.  F(  and  Ffy  are  the  forces  due  to  the  conservative  interactions,  while  F./J  and  F;f 

are  the  dissipative  and  random  forces,  respectively.  A  particle  temperature,  0t ,  is  defined  as 
Q.  =  duJdSj ,  where  .s/  is  defined  as  a  mesoscopic  entropy.  In  general,  this  microscopic  state  law, 
or  mesoparticle  equation  of  state,  6  =  0{u) ,  can  be  determined  from  molecular  simulations,  first- 
principle  calculations,  or  available  experimental  data. 

We  induced  steady  planar  shear  flow  by  means  of  the  Lees-Edwards  boundary  conditions  (6),  in 
which  the  simulation  box  and  its  images  centered  at  (x,  y)  =  (±L,  0),  (±2 L,  0),  ...,  are  taken  to  be 
stationary,  while  boxes  in  the  layer  above,  (x,  y)  =  (0,  L),  (±L,  L ),  (±2 L,  L ),  ...,  are  moving  at  a 
speed  tuL  in  the  positive  x  direction,  where  57  is  the  shear  rate.  Boxes  in  the  layer  below,  (x,  y)  = 
(0,  - L ),  (±L,  - L ),  (±2 L,  - L ),  ...,  move  at  a  speed  ujL  in  the  negative  x  direction.  A  system  under 
such  conditions  is  subjected  to  a  uniform  steady  shear  in  the  xy  plane. 

Next  we  parameterized  a  mesoparticle  potential  for  TNT  by  fitting  to  the  Mie-Griineisen 
equation  of  state  at  several  state  points.  Mesoparticles  interact  through  a  pairwise  additive  third- 
order  Rydberg  potential  given  as 


¥r 


j-e(l  +  d  +  ad})e  d 

I  0 


forfy  ^ 

for  rj>rcat 


(4) 


3 


.  As  a  first 


where,  r(/-  is  the  separation  distance  between  particle  i  and  j,  and  d  =  a 


,1 


approximation,  we  perfonned  fitting  by  specifying  ro=5.07  A,  o= 7.90,  a=0. 1 85,  and  rcut=  15.0  A, 

—20 

and  subsequently  detennining  £=1.612x10  J.  A  mesoparticle  was  chosen  to  represent  a  single 
TNT  molecule. 


From  statistical  thermodynamics,  the  internal  particle  energy  of  an  isolated  nonlinear  polyatomic 
molecule  containing  n  number  of  atoms  can  be  expressed  as  the  sum  of  the  translational, 
rotational,  vibrational,  and  electronic  contributions,  respectively,  so  that 
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where  &vib  .  is  the  characteristic  vibrational  temperature,  and  De  is  the  depth  of  the  ground 

electronic  state  potential  well  (7),  both  of  which  can  be  determined  from  quantum  mechanics 
calculations  or  taken  from  thennodynamic  tables  (5)  when  available,  and  where  kg  is  the 
Boltzmann  constant,  T  is  the  temperature.  As  a  first  approximation  for  this  study,  we  formulate 
equation  5  based  upon  the  classical  limit  (relatively  high  temperature)  (7),  so  that 

—  =  -  +  -  +  (3(21)-6)=60,  (6) 

kBT  2  2  V  ' 


since  each  TNT  molecule  contains  21  atoms  and  is  represented  by  a  single  mesoparticle.  By 
replacing  T  with  the  internal  particle  temperature  9i ,  we  arrive  at  the  equation  of  state  for 

particle  i 


9.  =■ 


60  kB 


(7) 


A  sub-grid  model  is  developed  to  account  for  the  dissipation  of  energy  and  allow  for  delayed 
reactions.  The  initiation  phenomenon  in  heterogeneous  energetic  materials  can  occur  when  the 
material  is  subjected  to  impulses  such  as  shock  waves.  Under  sufficiently  strong  shock 
conditions,  these  shock  waves  can  evolve  into  self-sustaining  detonation  waves.  This  process, 
shock-to-detonation  transition  (SDT),  is  fairly  well  understood  at  the  phenomenological  level 
and  is  based  on  the  theory  of  hot  spot  formation  (such  as  void  collapse,  visco-plastic  heating, 
shear  band,  frictional  heating,  etc.).  However  for  conditions  such  that  the  SDT  process  does  not 
occur,  the  initiation  of  HE  leading  to  explosion  can  take  place  if  the  material  has  been 
sufficiently  confined.  Such  analysis  was  shown  by  Frey  (P),  whose  study  indicated  that  both 
pressure  and  shear  rate  were  important  parameters  in  controlling  runaway  explosion.  Another 
study  showed  that  a  shear  banding  mechanism  could  provide  the  large  ignited  surface  area, 
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which  is  believed  to  be  necessary  to  explain  shock  initiation  (1).  The  long  delay  between  the 
time  at  impact  and  explosion  can  typically  be  hundreds  of  milliseconds;  whereas,  for  the  SDT 
process,  it  is  on  the  order  of  microseconds.  Such  dominant  physical  features  suggest  energy 
competing  processes  at  the  micro-scale  level.  Therefore,  at  hot  spot  locations,  e.g.,  shear 
surfaces,  the  possibility  of  initiation  hinges  on  a  balance  (or  lack  thereof)  between  energy 
producing  mechanisms  (visco-plastic  work,  shear  localization,  chemical  reaction,  etc.)  and  the 
rate  at  which  the  energy  is  transported  away  from  the  zone.  Energy  generated  from  shear 
localization  within  the  narrow  region  of  the  shear  band  width,  which  is  on  the  order  of 
micrometers,  can  be  much  higher  than  the  bulk  energy  in  the  typical  cell  size  of  continuum 
simulations,  which  is  on  the  order  of  millimeters.  The  large  temperature  gradient  within  a 
continuum  cell  supports  the  need  to  account  for  thermal  diffusion.  Figure  2  illustrates 
temperature  profiles,  including  energy  release  due  to  reaction,  as  time  progresses.  This  concept 
forms  the  basis  of  our  sub-grid  model. 

Within  each  sheared  continuum  cell,  we  assume  that  the  shear  band  surface  is  located  at  the  mid¬ 
plane  (figure  2).  Such  a  geometric  assumption  leads  to  a  simple  and  tractable  solution  at  the  sub¬ 
grid  level.  Each  sheared  cell  is  subdivided  into  intervals,  and  time  dependent  equations  of 
temperature  and  species  are  described  as  follow: 


pC  —  =  k  <—^-r  +  A/ Ini" 

dt  dx2 


dK 

p^t  =  pd 

dt 


dx 2 


-  m 


(8) 

(9) 


where,  X  is  the  extent  of  reaction  at  the  sub-grid  level,  heat  of  reaction  is  AH ,  diffusion 

coefficient  is  D ,  thermal  conductivity  is  k  ,  and  the  volumetric  mass  production  rate,  m" ,  is 
governed  by  an  Arrhenius-type  reaction: 


m =  Ap(l  -  Xsg  )exp^  1 j,  (10) 

where  A  is  the  frequency  and  Ta  is  the  activation  temperature.  The  species  diffusion  term  is 

assumed  negligible  in  the  current  study.  Solutions  for  the  above  equations  are  integrated  using  an 
implicit  method  with  the  initial  conditions: 
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where  the  shear  band  temperature,  Tmeso ,  is  obtained  from  mesoscale  simulations  of  a  shear 
banding  process,  Tcell  is  the  continuum  cell  temperature,  and  wsh  is  the  shear  band  width. 
Boundary  conditions  are  given  as 


dT_ 

d*  x=o 


(12) 


With  each  sheared  cell,  a  volumetric  averaged  solution  of  the  extent  of  reaction  at  the  sub-grid 
level  is  transferred  to  the  continuum  level  as  the  continuum  extent-of-reaction.  This  extent-of- 
reaction  is  used  in  the  equation  of  state  to  determine  the  amount  of  energy  release.  Such  an 
approach  ensures  consistency  in  our  internal  energy  calculation. 


3.  Results 


We  generated  the  mesoscale  look-up  table  for  a  range  of  densities  (1.65,  1.70,  1.80,  1.90,  and 
2.00  g/cm  )  and  three  shear  rates  (50,  75,  and  100  s  ),  calculating  both  pressure  and  temperature. 
Figure  3(a)  shows  a  sample  of  the  data  found  in  the  look-up  table.  Figure  3(b)  shows  a  sample 
energy  spectrum  determined  by  the  mesoscale  model  depicting  the  internal  temperature  spike 
due  to  shearing. 
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Figure  3.  (a)  Mesoscale  temperatures  calculated  as  a  function  of  pressure  and  shear  rate,  and  (b)  temporal 

variations  of  internal  (black)  and  kinetic  (red)  temperatures  after  the  system  has  been  sheared  in  the 
mesoscale  simulation. 


Using  TNT  parameters,  we  subjected  a  block  of  material  to  pure  shear  by  imposing  a  velocity  on 
the  lower  half  of  the  block  to  verify  the  newly  implemented  approach  in  CTH,  as  shown  in  figure 
4.  The  simulation  represents  a  hypothetical  configuration  to  simply  conduct  a  numerical  test  and 
demonstrate  the  depiction  of  a  reaction  that  previously  was  unattainable. 


Shearing 

surface 


(a) 


(b) 


(c) 


Figure  4.  Example  problem  to  verify  a  new  multi-scale  approach:  (a)  sliding  half  a  block  of  material,  (b)  extent  of 
reaction  at  2.5  ps,  and  (c)  extent  of  reaction  at  5  ps. 
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4.  Conclusions 


In  this  effort,  we  developed  a  multi-scale  approach  for  simulating  shear  initiated  reactions  that 
span  from  the  molecular  scale  to  the  mesoscale  (length  scale  of  material  heterogeneities)  to  the 
full  continuum  scale  (length  scale  of  the  weapon  system).  We  then  constructed  a  framework 
based  on  our  approach,  implemented  it  into  the  CTH  hydrocode  developed  by  Sandia  National 
Laboratory,  and  demonstrated  that  the  new  approach  allows  predictions  not  previously  possible 
for  energetic  materials  when  subjected  to  loads  that  result  in  shear  localizations. 

This  computational  tool  provides  a  novel  modeling  capability  that  opens  previously  unavailable 
avenues  by  bridging  the  gaps  between  multiple  scales.  This  capability  enables  improved 
predictions  towards  designing  armor  and  anti-armor  devices.  It  also  supports  the  development  of 
concepts  for  enhanced  survivability  and  lethality,  primarily  in  the  areas  of  insensitive  munitions 
and  reactive  armor,  and  of  novel  concepts  and  designs  using  RMs. 
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6.  Transitions 


Transition  of  the  multi-scale  model  to  Sandia  is  planned  and  will  be  executed  following  the 
completion  of  the  effort  in  fiscal  year  2009  (FY09). 
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List  of  Symbols,  Abbreviations,  and  Acronyms 


CTH 

DPDE 

FY09 

HE 

RM 

SDT 


Combined  Hydro  and  Radiation  Transport  Diffusion 

Dissipative  Particle  Dynamics  with  Energy  Conservation 

fiscal  year  2009 

high  explosive 

reactive  material 

shock-to-detonation  transition 
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Copies  Organization 

1  DEFENSE  TECHNICAL 
(PDF  INFORMATION  CTR 
only)  DTIC  OCA 

8725  JOHN  J  KINGMAN  RD 
STE  0944 

FORT  BELVOIR  VA  22060-6218 

1  CD  DIRECTOR 

US  ARMY  RESEARCH  LAB 
IMNE  ALC  HR 
2800  POWDER  MILL  RD 
ADELPHI  MD  20783-1197 

1  CD  DIRECTOR 

US  ARMY  RESEARCH  LAB 
AMSRD  ARL  Cl  OK  TL 
2800  POWDER  MILL  RD 
ADELPHI  MD  20783-1197 

1  CD  DIRECTOR 

US  ARMY  RESEARCH  LAB 
AMSRD  ARL  Cl  OK  PE 
2800  POWDER  MILL  RD 
ADELPHI  MD  20783-1197 

ABERDEEN  PROVING  GROUND 


6  HCs  DIR  USARL 

1  CD  AMSRD  ARL  Cl  OK  TP  (BLDG  4600)  (1  CD) 
AMSRD  ARL  WM  TC 
M  FERMAN-COKER  (5  HCs) 

AMSRD  ARL  WM  BD 
J  BRENNAN  (1  HC) 

TOTAL:  11  (1  ELEC,  6  HCs,  4  CDs) 
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